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We investigate few-boson systems in finite one-dimensional multi-well traps covering the full inter- 
action crossover from uncorrelated to fermionized particles. Our treatment of the ground state prop- 
erties is based on the numerically exact Multi-Configurational Time-Dependent Hartree method. For 
commensurate filling we trace the fingerprints of localisation, as the interaction strength increases, 
in several observables like reduced density matrices, fluctuations and momentum distribution. For 
filling factor larger than one we observe on-site repulsion effects in the densities and fragmentation 
of particles beyond the validity of the Bose-Hubbard model upon approaching the Tonks- Girardeau 
limit. The presence of an incommensurate fraction of particles induces incomplete localisation and 
spatial modulations of the density profiles, taking into account the finite size of the system. 

PACS numbers: 03.75.Hh,03.75.Lm 05.30.Jp 



I. INTRODUCTION 

Ultracold gases and Bose-Einstein condensates repre- 
sent a highly controllable and rich many-body system 
[l|, [2| , appropriate for the investigation of several quan- 
tum phenomena. The possibility to accurately design the 
external forces employing static or laser fields, along with 
the tuning of the interaction strength by visiting Fesh- 
bach resonances, allows us to prepare and explore weakly 
or strongly correlated quantum systems in almost arbi- 
trary potential landscapes [3]. 

In particular, the flexible experimental toolbox of opti- 
cal lattices enabled the investigation of the phase transi- 
tion from a coherent superfluid (SF) to a Mott-insulator 
(MI) state with localized particles in each site [^ , studied 
theoretically by the Bose-Hubbard model (BHM) [5|, |6| . 
Further highlights include quantum phases like Bose glass 
and Mott-shells, occurring for disordered, confined and 
incommensurate systems [Ml6 



Tuning appropriately the laser beams it is possible to 
effectively decrease the dimensionality of the system and 
explore intrinsic phenomena of quasi one-dimensional 
(ID) systems, where the so-called confinement-induced 
resonances resulting from the change of the transverse 
length scale [l7| can be utilized to alter the interaction 
strength. This way, it became experimentally feasible to 
investigate weak to strongly interacting systems and es- 
pecially the extreme case of the Tonks- Girardeau (TG) 
gas, a state of strongly correlated bosonic matter intrinsi- 
cally tied to ID physics [18|, [l9| where bosons that repel 
each other with infinitely strong forces can be mapped 
to non-interacting fermions [20]. Not only the special 
features of such a gas have attracted theoretical studies 
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but also the crossover to this limit from 3D to ID 



While the analytical treatment of the interaction 

crossover in ID from uncorrelated to TG-gas is possible 

for a homogeneous system with periodic [26] or hard- wall 

23 boundary conditions, the exact study of bosonic sys- 



tems exposed to specific trapping geometries is particular 
for a small number of particles. The harmonic trap [28[ 



311 and the double well 
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35I-I37J are paradigm sys- 
tems which unveil characteristic features of this crossover 
in confined geometries. Broadening and immergence of 
oscillations of the densities, as well as the evolution of 
the coherence properties have been examined in detail. 

ID optical lattices are very special and appealing sys- 
tems since they combine the physics in a lattice (phase 
transitions) with the special features of one dimension 
[sl, I22I, |38|-|42| . Localization, delocalisation as well as lat- 
tice imperfection effects have been studied for small en- 
sembles [43] showing interesting analogies with macro- 
scopic phases. In particular momentum distributions, 
pair correlations and energy spectra give insight to the 
MI-SF transition with increasing lattice depth or unveil 
the effect of an incommensurate filling, while a Bose-glass 
phase emerges when breaking the lattice symmetry. The 
effect of a higher filling factor and the subsequent on-site 
fragmentation of particles in periodic ID lattices has been 
explored in [4^] . The authors investigate a commensurate 
case with two particles per site and distinguish between a 
phase of MI with unperturbed Wannier functions, and a 
second transition to a state with two fragmented orbitals 
on each site. 

In this article we cover the different interaction regimes 
from the weakly correlated to the fermionization limit for 
a few-boson ensemble in a finite one-dimensional lattice. 
We go far beyond the BHM regime where the change of 
the interactions is equivalent to the change of the lattice 
depth. Since commensurability is a key issue concern- 
ing the crossover, we examine representative few-body 
setups for commensurate and incommensurate filling fac- 



tors. We demonstrate localization mechanisms by exam- 
ining densities, particle number fluctuations and coher- 
ence loss in momentum distributions in the case of a unit 
filling factor. For higher commensurate filling we observe 
on-site repulsion effects such as broadening and oscilla- 
tory patterns in the densities which go beyond the Bose- 
Hubbard regime. Incommensurable setups on the other 
hand, show the possibility for spatial variations of observ- 
ables, as a result of the correlations and the finite size of 
the system, in addition to partial delocalisation due to 
the incommensurate fraction of particles. The method 
which we use here is the Multi-Configurational Time- 
dependent Hartree (MCTDH) [45i . |46| which has proven 
to be efficient in treating bosonic systems both for static 
properties and dynamics (see Appendix) [36, 37, 47|. 

The paper is organised as follows: In Sec. II we explain 
the setup, the analytically describable limits and the ap- 
proximative methods. The results are presented in Sec. 
Ill for commensurate and Sec. IV for incommensurate 
filling. We summarise our findings and give an outlook 
in Sec. V. A short description of the numerical method 
(MCTDH) is included in the Appendix. 



II. SETUP AND THEORETICAL 
BACKGROUND 

A. Model Hamiltonian 



cases of commensurability. We use a sufficiently large 
lattice depth {Vq = 7.0 — 20.0 which is of the order 
of 4 — 20 Er depending on hi) such that at least two 
single particle bands lie energetically below the contin- 
uum for reasons we explain later on. The way we render 
our system finite via imposing hard- wall boundaries does 
not restrict the generality of our treatment and results. 
In particular, independently of the specific experimental 
implementation (using multi-colour lattices and/or har- 
monic confinement), finite systems exhibit intrinsic fea- 
tures which differ from infinite lattices (periodic bound- 
ary conditions); the confined traps result always in spa- 
tial inhomogeneities of the densities which are absent 
in the periodic case. These inhomogeneities can be in- 
creased or manipulated if there is a harmonic confinement 
or a disordered surface. Although we are focusing here 
on the case of equal on site energies, we observe a rich 
behaviour that captures the main effects of a finite con- 
fined system. All our numerically exact calculations are 
performed by the MCTDH method (see appendix for a 
description). The following subsections referring to the 
extreme cases of no and infinite interactions, as well as 
the Bose-Hubbard model for weak interactions, are valid 
only in certain regimes which we address in the Sees. HI, 
IV to explain the observations. Covering the complete in- 
teraction crossover from ^ = to ^ ^ oo by MCTDH we 
are able to show the regimes of validity but also examine 
effects going beyond these models. 



The effective ID Hamiltonian we consider reads: 



N 
H = ^hi^^ Vint{Xi - Xj) 
i i<j 



(1) 



where the two-body interaction potential is delta-like 
Vint{xi-Xj) = giD^Xi-Xj). The coupling strength qid 
depends on the scattering length ao and the oscillator 

length ax = 



h 

Muj 



which characterises the transverse 



confinement: giD = 1^(1 - ^73^)"^ 113, ll^- The 

2 
one-body part of the Hamiltonian hi = -^ ^V{xi) con- 
tains the ID- lattice potential V{xi) = Vq ^[^^{tvXi)^ char- 
acterised by the depth Vq and periodicity d (distance be- 
tween two successive minima) setting k = ir/d. In order 
to restrict the infinite potential V{xi) = Vq sin^(/^x^) to a 
finite number of sites W and a length L, we impose hard- 
wall boundary conditions on appropriate position. In our 
calculations, for numerical convenience (keeping a stan- 
dard grid), the length unit is ay = L/10 and the energy 
unit is chosen as fi^/Man setting also h = M = 1. The 

rescaled scattering strength is thus g = ^" ^2^^ . In units 
of the recoil energy Er - 



2 ..2 



2M 



it reads g' = j^ = 



2d 



a\\ TV' 
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All the parameters are considered to be tunable almost 
at will in corresponding experiments. For the interac- 
tion strength g we take representative values to cover 
the complete crossover to fermionization in all different 



B. Single Particle states 

In order to understand the limit cases of non- 
and infinitely- interacting particles (see next subsec- 
tion), a discussion of the single particle states is 
necessary. Analytical expressions for the delocalised 
single-particle states, i.e. Bloch states, are avail- 
able for periodic boundary conditions. For finite lat- 
tices, we use the tight-binding approximation assum- 
ing only a nearest-neighbour tunneling coupling term 
J (X — J Ws{x)hiWs-\-i{x)dx between the sites s and 5 + I, 
where Ws{x) are the on-site localised Wannier states. 
Within this approximation, valid for a relatively deep 
potential, we express the Bloch states in terms of Wan- 
nier functions. The single-particle Hamiltonian written 
in a matrix form using the localised states basis is: 



hi = 







-J 

62 -J 
-J 63 



\ 



:/ 



where e^ {s = 1,2,..,VK) are the on-site energies, which 
in our case are equal (ei = ... = ew = e). 

The hard- wall boundary conditions imply that there is 
no tunnel coupling between the first and the last lattice 
site as opposed to the periodic ones, where there is a 
single coupling for all sites. The resulting eigenvalues 




Figure 1: (Colour online) Sketch of the finite three- well lattice 
and the corresponding single particle states for the first two 
bands. 



are: Eq-i = e — 2Jcos( 
eigenfunctions read: 
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(2) 



These lowest-band single-particle eigenstates are 
sketched in Fig. [1] for the triple well using the Gaussian 
approximation for the Wannier functions around the 

center of each well Xg {x\ws) = (7r(i^)~"'^/^e~^^~^^^ /^^ : 
(fo) = ^(ki) + V2\w2) + 1^3)), \^i) = 75(ki) - ks)), 
^2) = ^(I'^i) — a/2 I '^2) + \^s))i with corresponding 

energies Eq = e — a/2 J, Ei = e, £^2 = e + V^J. Note that 
for the ground state the middle well is occupied with a 
larger amplitude compared to the two outer ones. For 
the states of the excited bands the harmonic oscillator 
orbitals of higher order can serve as localised functions 
to a rather good approximation. The computation of 
the eigenstates, done by numerical diagonalisation of 
the model Hamiltonian hi is in good agreement with 
this simple model h^. In a recent paper [13], the single 
particle states of the lowest band for the case of a 
lattice with a superimposed parabolic trap were derived 
analytically within the tight-binding approximation. 



C. Bose- Fermi map 

The non-interacting and the infinitely strongly inter- 
acting case of impenetrable bosons can be both addressed 
analytically, given that we can solve the corresponding 
single-particle problem. In the former case, we just have 
one orbital for the ground state (the lowest eigenstate 
of the single-particle Hamiltonian) where all bosons re- 
side "^0 = (J)q . In the Tonks- Girardeau limit, the 
Bose-Fermi map [20] establishes the following isomor- 
phy between the bosonic ^tg and the non-interacting 
fermionic wave function ^fer* ^tg(Q) = ^(Q)^fer(Q) 
where A{Q) = ni<jSgn(xi — Xj) and Q = (xi, ...,XAr)"^. 
In particular, the infinitely interacting bosonic ground 
state is simply the absolute value of the fermionic one. 



The local densities are appropriate tools for analysing 
a many body state |^): the one body density p{x) = 
{x\pi\x)^ diagonal kernel of the one-body density op- 
erator pi = tr2,...,Ar|^)(^|, and the two body density 
^2(^1,^2) = {xiX2\p2\xiX2) diagonal kernel of the two- 
body density operator p2 = tr3^...^Ar|^)(^|. In the TG 
limit they are equal to the fermionic ones: 



N-l 



p^^'i^) = Yl i^-(^) 



a=0 



T.O<a<a^<N-l \^a{xi)(Pa'{x2) - (Pa{x2)^a' {xi)\^ . 

Accordingly the ground state energy in the fermioniza- 
tion limit is ^tg(A^) = ^a=o ^oc- 



D. Bose-Hubbard model and beyond 

Bosons in optical lattices have bee mostly studied in 
the literature within the Bose-Hubbard model [3], em- 
ploying the tight-binding and lowest band approxima- 
tion. Apart from the coupling term J and the on-site 
energies e^ introduced in hi^ there is additionally on-site 
interaction with strength U which reads for a delta con- 
tact potential U = g j dx\ws{xy^' 
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= -J Y^ alas' ^izY 
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where < 5,5^ > indicates the sum over nearest neigh- 
bours. Both parameters J and U can be tuned by the 
lattice constants Vq and d. Increasing for example the 
laser intensity, the lattice becomes deeper, and thus the 
ratio U/J increases substantially leading to the MI phase 
in the commensurate case: the particles save interaction 
energy by being localised in different wells and the tun- 
neling to neighbouring sites is strongly suppressed. The 
SF phase on the other hand, is characterised by phase 
coherence and delocalisation of the particles. Quantum 
phase transitions at zero temperature are triggered by 
quantum fluctuations, and in this sense, we can exam- 
ine in our few-body ensembles signatures of the MI and 
SF phase. We use as a measure of localisation, the local 
(on-site) particle number fluctuations: 



ANg^ 



{ug') - {ugf = N[p2AN 



l) + pi,(l-ArpiJ] 

(3) 
where pi^ = J^dxp{x) and p2, = J^ Jgdxidx2p(xi,X2) 
are the one- and two- body densities respectively, in- 
tegrated over the lattice site s. The SF-MI transition 
is accompanied by decreasing and finally vanishing fluc- 
tuations ANg as U/J increases (corresponding in some 
sense to vanishing local compressibility in the MI phase 



[3). Another sign of the transition is the loss of coher- 
ence due to the locahsation of the particles in individ- 
ual sites. This can be observed in the momentum dis- 
tribution, where the visibility of the interference peaks 
in the coherent SF phase is reduced, ending up with a 
smoothened incoherent profile in the MI regime [4']. 

For filling factors higher than one, i^ = -^ > 1, it 
is necessary to go beyond the simple BHM, which as- 
sumes unperturbed Wannier orbitals and is restricted 
to the lowest band, one has to include higher band ef- 
fects [43| to examine the fermionization and generally 
strong correlation effects. The Bose- Fermi map already 
indicates that the lowest band levels are not sufficient 
to accommodate a number of fermionized bosons larger 
than the number of wells. Let us briefly mention effective 
models relevant for the discussion of filling factors higher 
than one. A modulation of the Wannier functions to take 
into account on-site interaction effects has been proposed 
48 , 49] or splitting into two orbitals into the same well 



44| . Another recent few-body study suggests to optimise 



the BHM parameters such that they agree with the ex- 
act results specifically in the strongly interacting regime 



50|. Besides this, the concept of extended fermionization 



valid in the extreme limit of BHM U/J^oo under the 
assumption that the particles occupy different layers of 
MI and SF character [51] gives a valuable picture for situ- 
ations where an incommensurate fraction of particles sits 
on a commensurate localised background. We emphasise 
that the following results are obtained by the numerically 
exact MCTDH method and we only refer to other models 
for explanation and comparison. 



III. COMMENSURATE FILLING 

We explore here firstly the commensurate filling case 
showing the fingerprints of localization for filling factors 
1^=^ = 1 and i^ = 2, the latter being particularly 
interesting for the study of on-site interaction effects with 
two particles per site. 



A. Filling factor z/ = 1 

a. The spatial distribution of the particles For the 
non-interacting ground state, the density of particles is 
larger for the middle sites and decreases as we go to the 
outer ones, illustrated here for 6 wells and 6 particles in 
Fig. [2fa) for ^ = 0. This occurs for setups with hard- 
wall boundary conditions but is generally characteristic 
of finite lattices (independent of the number of sites and 
particles), since the kinetic energy term renders the mid- 
dle wells energetically more favourable. We remark that 
the particles reside all in the same single particle ground 
state which has indeed the maximal density in the mid- 
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Figure 2: (Colour online) (a) One body-density p{x) for 6 
wells and 6 particles. Shown are 4 different values of g: non 
interacting (g — 0.0), weakly interacting {g — 0.2,0.6), Mott 
insulator- fermionization limit [g — 10.0). (b) Particle num- 
ber fluctuations as g increases: red, green and blue lines cor- 
respond to the left 3 wells of the potential s — 1, 2, 3 counting 
from the outermost one. (inset) On-site population as g in- 
creases. 



serve a gradual redistribution of the density which leads 
to equal population of all sites. Let us explain this ob- 
servation in terms of the BHM which predicts a sim- 
ple localisation process where each particle sits in one 
well to save interaction energy. In the Fock state rep- 
resentation \Ni^N2^ '"^Nw)^ where each basis vector is 
parametrised by the occupation numbers for each site, 
with increasing /7/J, the vector which has no double 
occupation 1 1,1,1,...) becomes the lowest eigenstate of 
the BHM Hamiltonian. As soon as the particles localise 
one per well, the increase of the interaction does not af- 
fect them anymore. In this fermionization limit their 
energy actually saturates to the fermionic value. Here 
J is fixed (since there is a fixed depth Vq = 12.0 corre- 
sponding to 6.2Er with d = 1.6) and U varies with g. 
In this special case u = 1 and in general for v < 1, the 
Bose- Fermi map does not hold only for the abstract state 



^TG < — > ^fer but also for the Wannier state l8|, |38|, |39 



As the interaction strength increases (Fig.[2fa)), we ob- 



In other words, the U/J -^ oo limit of the BHM where 
each particles sits in one well |1, 1, 1, ...), is here equiva- 
lent to the fermionization limit, where each particle oc- 
cupies a single-particle level. Thus the BHM is valid in 
the case of z^ = 1 for the whole range of interactions. 
In Fig. [2](b) (inset) we plot the population of each well, 
Npi^ , as a function of the interaction strength. While for 
the outer (5 = 1) and middle {s = 3) wells the population 
evolves monotonically to the final value, the inter medi- 



ate wells (see 5 = 2) act as a 'carrier' of particles and 
thus their population may exceed for some intermediate 
interaction strength {g ^ 0.1) their final one. 

b. Diminishing fluctuations as a sign of localization 
The uniformisation of the population with increasing 
interaction strength is a pre-signature for the local- 
ization. A more accurate measure indicative for this 
mechanism, is the particle number fluctuations (Eq. 3) 
which decrease substantially as the interaction increases 
(Fig. Efb)) . Of course, a lower J, corresponding to a 
deeper lattice, enhances the localisation process. For a 
very deep lattice U/J -^ oo we expect them to vanish 
completely, while in our case they saturate to a small 
value. The rather shallow depth of our lattice permits an 
occupation of the interwell space with a non-zero overlap 
between particles sitting in neighbouring sites. Includ- 
ing this area into the integrations of the one- and two- 
body density (Eq. 3) we end up with a small contribu- 
tion to the fluctuations even in the strongly interacting 
limit. Note also that the middle wells {s = 3), which 
can 'lose' population with respect to both sides, keep on 
having larger fluctuations than the outer ones {s = 1). 

In Fig. El^b) we observe that the fluctuations converge 
to a constant value for g ^ 3.0. This convergence hap- 
pens for values of g signiflcantly larger than those where 
the uniform distribution of the population is achieved 
{g = 0.6 Fig. [2fb) inset). Technically speaking, particle- 
hole excitations like |2,0,1,...), |0,2,1,...) contribute for 
g ^ 0.6 resulting in an equal site distribution but without 
forming a localised state |1, !,..,!). Perfect localisation 
occurs only when the latter vector is the eigenstate of the 
system suppressing all other contributions. 

Let us point out here the effect of the hard wall in 
comparison with periodic boundary conditions. In our 
conflned system for weak interactions, the Fock states 
with less occupation in the outer wells outweigh those 
with less occupation in the center and therefore we ob- 
serve an imbalance of the sites population. In the case 
of periodic boundary conditions even the non-interacting 
state possesses a uniform population of the sites due to 
symmetry. Hence, what really happens in the latter case 
within the transition from weak to strong interactions, 
is only a reduction of fluctuations by going to the state 
|1,1,1,...), and not a redistribution of the site popula- 
tions as in our case. 

c. Two-body correlations In the two-body density 
p2{xi^X2) (Fig.[3j), starting from g = and for increasing 
interactions {g = 0.2), the diagonal peaks indicating dou- 
ble occupation (especially the middle) fade out, while the 
off-diagonal ones are amplifled. Note that in the regime 
where equal distribution over the sites has been achieved 
{g = 0.6), the diagonal of the two body density has not 
yet been fully emptied, a sign that some double occupa- 
tion is still present and thus a true localised state with 
unit fllling is not yet reached. As ^ ^ oo the repulsive 
forces totally prohibit double occupation and lead to a 
complete depletion of the diagonal {g = 10). As a result, 
if one particle sits in one well, any second one distributes 
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Figure 3: (Colour online) Two body-density p2{xi,X2) for 6 
wells and 6 particles. Shown are 4 values of the interaction 
strength (a) g = 0.0, (b) g = 0.2, (c) g = 0.6, (d) g = 10.0. 



itself equally over the other wells, but has zero probabil- 
ity to be in the same well. 

d. Non-local correlations We will now cast light on 
the system from the perspective of the non-local prop- 
erties speciflcally the off-diagonal kernel of the one-body 
density matrix pi{x^x') = {x\pi\x') and the momentum 
distribution. While the local properties in the TG-limit 
have exactly the fermionic proflle, the bosonic permu- 
tation symmetry plays a signiflcant role for the coher- 
ence properties. The off-diagonal behaviour of pi{x^x') 



as 



x—x'\ -^ oo is a measure of coherence [52|- it indicates 
non- vanishing off-diagonal long range order (ODLRO) 
in inflnite homogeneous systems. In our flnite setups 
though there can be no true ODLRO, so the term co- 
herence refers here to the off-diagonal parts of piix^x') 
showing short and long range one-particle correlations. 

For zero interactions (see Fig. H] ^ = 0.0) we observe 
that the off-diagonal spots fade out for an increasing dis- 
tance from the center as expected for a conflned system. 
For weak interactions {g = 0.2), the off-diagonal con- 
tributions -especially the remote ones- become more pro- 
nounced, along with the outer wells on the diagonal. This 
is a common feature for flnite setups, being accompanied 
by the initial density redistribution. The particles, re- 
distribute all over the lattice such that they reduce the 
interaction energy, and as long as they stay in one orbital 
(see explanation later on), increase the correlations all 
over the space. As the interaction strength increases fur- 
ther {g = 0.6), the localisation on discrete sites gradually 
destroys the ODLRO . However, some short range coher- 
ence persists as the distribution of the particles becomes 
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Figure 4: (Colour online) Off diagonal one-body density ma- 
trix pi{x, x') for 6 wells and 6 particles. Shown are 4 values of 
the interaction strength (a) g — 0.0, (b) g — 0.2, (c) g — 0.6, 
(d) g = 10.0. 



uniform {g = 0.6) and dies out for a stronger value of the 
interaction strength {g = 10.0). This matches with our 
previous observation, that the diminishing of the fluctua- 
tions and coherence in the MI phase does not necessarily 
coincide with the appearance of a uniform population. 

Momentum distribution The off-diagonal part of 



e. 



the one-body density matrix is not an observable but it is 
indirectly accessible via time-of- flight measurements J3|, 
which yield the momentum distribution p{k): 



p{k) = 27r{k\pi\k) = dx dx'e '^"^^ ^'^pi{x,x'). 



(4) 



In Fig.[5fa) we observe that p{k) exposes a rich pattern 
for ^ = with Bragg peaks near the reciprocal lattice vec- 
tor (/c = a* = ^ (i = 1.6) which is gradually smeared 
out as the interaction strength increases. The central 
peak corresponding to the ODLRO is high but gets even 
slightly higher for small interactions {g = 0.05), match- 
ing with our observation in pi(x, x') that the remote off- 
diagonal humps increase. In this SF coherent regime 
{g = 0.0,0.05,0.2) there occur also minor dips at the 
points km = (m/W)a*, due to the suppression of stand- 
ing waves with odd parity and wavelengths A^ = ^^ 



m 



{m = 1,2,...) resulting from the confinement, since the 
ground state possesses an even parity [43]. For ^ = 0.6 a 
value corresponding to uniform distribution in the den- 
sity, the central peak is substantially lowered, which is a 
typical sign of coherence loss because of localisation, but 
the lattice geometry fingerprints (side peaks) in the mo- 
mentum profile do not vanish completely implying again 



imperfect localisation. As the interaction increases fur- 
ther {g = 10.0), p{k) goes to a smoothened gaussian-like 
profile with complete destruction of the superfluid inter- 
"erence pattern and the visibility of the associated peaks 
^, |43|. We end up with a 'MI- type' incoherent state. 



/. Fragmentation analysis via natural orbitals The 
complete information in the one-body level is given by the 
spectral decomposition of pi = ^i=Qni\4'i){4^i\^ where 
the relative populations ni serve as a measure of frag- 
mentation into effective single particle states (pi (natural 



orbitals). In [53J, a criterion for the a non-fragmentated 
condensed state was introduced exactly by demanding 
the highest such occupation no to be very close to one. 



For zero interactions, all the particles reside in the low- 
est (ground state) natural orbital (see Fig. [51(b) g = 0.0). 
Increasing the interaction they gradually fragment into 
the first N = W = 6 orbitals. Thus, in this case 
{u = 1) the effective description through the lowest-band 
single-particle states holds. The population of each or- 
bital in the fermionization limit is not exactly 1/A^ as we 
would naively expect from a mapping to non-interacting 
single-particle states according to Girardeau's theorem 
(for fermions ni = 1/7V for I = O^N — 1). This is because 
the natural orbitals are effectively modulated single par- 
ticle states originating from the spectral decomposition 
of pi, and this modulation accounts for interaction ef- 
fects. The lowest orbital (/)o, for example (see Fig. [5fc)), 
is broadened for weak interactions {g = 0.2), follow- 
ing the evolution of the one-body density. The latter 
fact validates the use of the Gross-Pitaevkskii mean- field 
treatment for weak interactions where all the particles 
are assumed to reside in one variationally modulated or- 
bital. The modulation of this dominant orbital is respon- 
sible for the initial extension of the off-diagonal range in 
pi(x, x') and consequently leads to an increase of the cen- 
tral peak in p{k) (see Fig. [51(a)). The fragmentation ob- 
served for higher interactions (see Fig. [S^b)), is beyond 
the regime of validity of the Gross-Pitaevskii equation, 
and due to the admixing of higher orbitals coherence is 
destroyed. This may be compared with the cases of the 
double well and the harmonic trap [37|: in the double 
well, the interaction immediately bridges the gap within 
the lowest-band doublet and thus destroys the coherence. 
Note that the ground state orbital in case of the double- 
well has already an equal distribution for the two wells 
and thus there is no dramatic flattening of the density due 
to the interactions as in the case of more wells examined 
here. In the harmonic trap there is an initial extension of 
the long-range order exactly for the same reason as here, 
i.e., the ground state orbital is broadened and holds the 
main population (the gap between the ground and the 
excited state is relatively large here). For very strong in- 
teractions (Fig. [SJc) g = 10.0), the profile of the orbitals 
tends to return to the non-interacting one indicating the 
validity of the Bose- Fermi map in the TG-limit. 
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Figure 5: (Colour online) (a) Momentum distribution for 6 
wells and 6 particles. Shown are 5 different values of g: non 
interacting (g — 0.0), weakly interacting (g = 0.05,0.2,0.6), 
Mott insulator- fermionization limit (g — 10.0). (b) Popula- 
tion of the natural orbitals as a function of the interaction 
strength, (c) Profile of lowest natural orbital for several val- 
ues of the interaction. 
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Figure 6: (Colour online) (a) p{x) for 3 wells and 6 particles. 
Shown are 5 different values of g\ non interacting {g — 0.0), 
weakly interacting [g — 0.02, g — 0.2), and on-site fermion- 
ization crossover [g — 5.0, g — 20.0). (b) Particle number 
fluctuations as g increases for the left (s = 1) and the middle 
(s = 2) well, (inset) On-site populations. 



B. Filling factor 1^ — 2 

In the case oi v = 2, here examined for 6 particles in 
3 wells, we expect on-site interaction effects to be im- 
portant because of the higher number density of parti- 
cles. The one-body density redistributes resulting to an 
equal population for all wells with increasing repulsion 
(see Fig. [6fa) g = 0.2). In terms of BHM, we have a 
formation of a 'Mott state' of 2 particles per site resid- 
ing in unaltered Wannier orbitals as the vector |2, 2, 2, ...) 
becomes eigenstate of the Hamiltonian for U/J -^ oo. 

After this localisation into pairs is achieved, on site 
interaction effects become apparent in the one-body den- 
sity beyond BHM. The well-studied and analytically solv- 
able [28] fermionization pathway of 2 particles in an in- 
dividual harmonic trap arises on each site. In partic- 
ular, there is a broadening of the one-body density in 
each well resulting from the increase of the onsite repul- 
sion (Fig. [6ja) g = 5.0). For even stronger interactions 
{g = 20.0), we observe the formation of two density max- 
ima per site . These patterns arising in the few body 
setups for the strongly interacting regime have been ex- 
amined for the harmonic trap [29, 31, 37], the double well 



34l-l37lJ54l and ID lattices with periodic boundary condi- 



tions 4^. In the latter case the authors distinguish two 



phases: first the localisation into pairs (BHM regime), 
and second fragmentation of each pair into two orbitals 
in the same well, which is consistent with our exact re- 



sults. We underline that the on site interaction effects, 
go beyond the validity of the BHM, because higher band 
contributions need to be taken into account. In terms of 
the Bose- Fermi map, for u = 2 all the levels of the first 
two bands are occupied in the fermionization limit. The 
upper band (see Fig. [T]) involves excited functions with 
one node per site. The combination of these functions 
with the lowest band ones, gives the observed wiggled 
fermionization profile with a local minimum in the mid- 
dle of each well. Modulation of the Wannier functions 
has also been proposed to include on-site two-particle in- 
teraction effects J48|. The particle number fluctuations 
(Fig.[6fb)) together with the populations (Fig.[6fb) inset) 
illuminate the whole process: first the particles distribute 
homogeneously {g = 0.0 — 0.1), then localise and the fluc- 
tuations tend to a low value (^ ^ 0.1 — 1.0) as expected 
for commensurate filling in the MI phase, and finally each 
pair of particles exhibits an individual on-site two-body 
crossover as described above. 

In the two-body density P2(^i, ^2) the diagonal contri- 
bution is reduced with increasing value of g (see Fig. [71 
^ = 0.2). In fact the maxima on the diagonal acquire half 
the population of any off-diagonal one, which matches 
with the pair localisation process leading to an equal pop- 
ulation of all sites. For stronger interaction (see Fig. [71 
^ = 5.0) the formation of a correlation hole at xi = X2 in 
the diagonal occurs. This is an inherent two-body effect 
of the on-site fermionization process which is smoothened 
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the deeper lattice that we use here compared to the case 
z^ = 1 in the previous subsection {Vq = 12.0 = Q.^En for 
z/ = 1 and Fo = 7.0 = 15.4^;^^ for z/ = 2 with d = 3.3). 
The fragmentation process is enhanced and thus the co- 
herence is directly destroyed by admixing higher orbit als. 
Indeed, in a heuristic single particle picture, the energy 
levels within one band come closer for a deeper lattice and 
thus the gap is bridged easily by the interactions. The 
instant approach of the populations of the contributing 
natural orbitals of the lowest band (see Fig. [8] (b) 0,1,2) 
accounts for this fact. The population of the natural 
orbitals reflects the band structure: orbitals 0-2 of the 
effective first band and 3-5 of the second band (being in- 
distinguishable in Fig. [8](b)). The former orbitals (0-2) 
have the dominant contribution and are modulated, fol- 
lowing quite well the evolution of the one-body density, 
including the on site interaction effects (see Fig. [8](d)). 



C. General remarks and energy properties 



Figure 7: (Colour online) p2{xi, X2) for 3 wells and 6 particles 
for (a) g = 0.0, (b) g = 0.2, (c) g = 5.0, (d) g = 20.0. 



out in the integrated one-body density. The 'incomplete' 
correlation hole is a significant characteristic of z^ > 1 
filling factors in general, as the particles have no chance 
to be in completely different wells and thus are obliged 
to minimise their overlap on the same site (diagonal). 
The correlation hole is a very prominent two-body effect 
and occurs even for interactions interaction strengths (at 
g ~ 0.6) where the on-site broadening and maxima in 
the one-body density are not yet pronounced. Neverthe- 
less, the formation of the correlation hole begins after 
the localisation of two particles per well is established. 
Also visible in Fig. [71 is a broadening (for g = 5.0) and a 
fragmented pattern (for g = 20.0), which appear in the 
off- diagonal. 

For pi(x, x') the main effects are similar to the case of 
V = 1. We comment that for the fermionization limit, 
pi{x^x') (Fig. [8] (a), g = 20.0) reflects the isolated max- 
ima of its diagonal part pi{x,x) = p(x), while the off- 
diagonal shows slight short-range correlations because 
of the broadening of the on-site functions. In the mo- 
mentum distribution (Fig. [8] (c)), the interference peaks 
(Bragg and central one) become lower already for the on- 
set of interactions, ending up in a complete smoothening 
for strong interactions g = 5.0, 20.0. The TG profile of 
the density {g = 20.0) differs from the localised state of 
the BHM {g = 2.0) in particular by the fact that the 
high momentum tails are more pronounced in the former 
case. Note that from ^ = 2.0 to ^ = 5.0 there is an in- 
crease of the k = peak which can be attributed to the 
on-site broadening of the density (similar to the case of 
the harmonic trap, see [36])- The immediate lowering of 
the central peak for small interactions is attributed to 



Let us generalise our findings before we move to the 
case of an incommensurate filling. The equal distribu- 
tion of the density onto the sites and the subsequent 
loss of fluctuations shown here, are also predicted by the 
BHM. These processes, which happen within the lowest 
band (and with Wannier states almost unaltered), are 
enhanced for a deeper lattice. On-site interaction effects 
occurring only for v > 1^ i.e, density broadening and for- 
mation of maxima as well as correlation hole and frag- 
mented patterns in the two-body density, only show up 
for strong interactions (beyond BHM) when a substantial 
population of particles is well localised in one site. This 
fermionization crossover also applies to any integer filling 
factor, with the on-site phenomena involving v = 2, 3, 4... 
particles per well. The off-diagonal one-body correla- 
tions increase for increasing but weak interactions with 
the particles dominantly occupying the first orbital which 
flattens initially. For stronger interactions the fragmen- 
tation of particles destroys the coherence and the high 
visibility of the peaks in the momentum distribution is 
washed out to a smooth profile. For two particles per 
site the high momentum tails show differences between 
the 'BHM insulator' and the Tonks limit. 

A general comment about the behaviour of the energy 
is in order here. The ground state energy increases with g 
and saturates for ^ ^ 00 to the corresponding fermionic 
one, as we have laid out in Sec. II C. One particular 
aspect of this crossover is the response of the energy when 
we switch on the interactions. Close to ^ = the slope 
for the energy is approximately 



dE 
dg 



N{N - 1) 



g=0 



(po{x)\^dx 



(5) 



For hard-wall boundary conditions we know from Eq. (2) 
that with increasing the number of wells (the size of the 
system), the population in the center of the lattice also 
increases if we keep the filling factor constant. Thus one 
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Figure 8: (Colour online) (a) pi(x,x) for 3 wells and 6 
particles in the fermionization limit g = 20.0. (b) Popula- 
tion of the natural orbitals as a function of the interaction 
strength, (c) Momentum distribution: shown are 4 differ- 
ent values of g: non interacting (g — 0.0), weakly inter- 
acting {g — 0.02,0.2,0.5), on site fermionization crossover 
{g — 5.0,20.0). (d) Profile of the lowest natural orbital for 
several values of the interaction. 



would expect that the repulsive interaction affects larger 
systems more strongly due to the existence of areas with 
higher density. On the other hand, the integral of the 
single-particle ground state wave function to the fourth 
power representing the interaction term in the Eq. 5 
above, is lower for a potential with more wells as the 
delocalisation of the state increases. The above effects 
cancel out resulting in an almost size-independent evo- 
lution of the ground state energy per particle as long as 
N = W [see Fig. [Ifa)]. Of course, for increasing filling 
factors interaction effects manifest more strongly in the 
evolution of the energy as the density plays the dominant 
role. 

A second aspect is the energy gap between the ground 
and the first excited state. Let us point out that even 
without interactions, we have here an energy gap between 
the two states of the order of J [see Fig. [9fb)]; a contin- 



uous band structure which results in a gapless spectrum 
in the SF regime arises only in the limit of an infinite 
lattice. Moreover, since the Bose-Fermi map holds also 
for excited states, we can compute exactly the gap for 
^ ^ oo, which is of course the interband gap of the sin- 
gle particle spectrum. In the BHM regime (see Fig. [9l^b) 
inset), the gap is of the order of [/, but in general, as 
considered here, it is evolving continuously and occurs 
not as a sudden transition as it does in the macroscopic 
case. 

The last but important comment refers to the role of 
the kinetic energy, which is related to the parameter J 
and this, in turn, to the lattice depth. We have already 
underlined its impact on the enhancement of localisa- 
tion and fragmentation for a lower J. We performed 
all the calculations assuming a sufficiently deep lattice, 
such that at least two single particle bands lie below the 
energy maxima of the barriers (see Fig. [1]). This choice 
validates our argumentation in terms of the tight-binding 
approximation and also it ensures validity of the BHM for 
small interactions. In the other limit of a shallow lattice, 
a hydrodynamic approach in the framework of the sine- 
Gordon model has been employed [8]. They show that, 
for very strong interactions in ID, an arbitrarily small 
perturbative lattice potential is enough to result in an 
insulating phase for commensurate filling. However, the 
latter discussion is based on the thermodynamic limit. 
For our system, if we have particles delocalised above the 
barriers, then the on-site few-body effects smoothen out, 
and there is an interaction induced broadening leading 
to a filling of the space between the wells. In general as 
the total energy increases the particles come energetically 
closer to the continuum, where the barriers thin out and 
thus there is an enhanced penetration into the barrier 
(flattened on-site functions) which effectively results in 
larger fluctuations. A hint towards this effect is already 
evident in the case of i^ = 2 in the one body density 
[Fig. [6] (a)] where for strong interactions g = 5.0, 20.0 
we see a slightly higher density in the inter- well space. 
Accordingly the fluctuations [Fig. [6] (b)] show a slope to 
slightly higher values . 



IV. INCOMMENSURATE FILLING 



Incommensurate filling is more susceptible to the ex- 
act number of particles compared to the number of wells, 
with one main feature: there is always a delocalised frac- 
tion of particles. We consider essentially two cases of u 
non-integer: u < 1 where on-site interaction effects do 
not manifest due to low population and v > 1 which 
for strong interactions can be interpreted as a fraction 
N mod W of extra delocalised particles sitti ng o n a com- 
mensurate background of localised particles [5l| 
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Figure 9: (Colour online) (a) Ground state energies per parti- 
cle as the interaction increases for commensurate filling i/ = 1 
for the cases 4 particles in 4 wells and 6 particles in 6 wells. 
(b) Two first eigenstates of the spectrum of 3 particles in 3 
wells as the interaction strength increases, (inset) weak inter- 
actions 



A. Filling factor u < 1 

The main concern here is how the particles distribute 
over the lattice as the repulsion increases. In the weak 
interaction regime, the repulsive forces drive the parti- 
cles away from the highly populated center of the po- 
tential, (see for example 5 particles in 7 wells in Fig. [10 
(a)). In this case the one-particle density does not tend 
to an equal site occupation, but stays even for strong 
interactions asymmetric {g = 3.0). The exact number 
of particles N and wells W determines the fermionized 
distribution of the density: 



w N 



VV 1\ / 

p(^)Fermi ^C ^ ^ | siu ( 

s q ^ 



sqiT 



W^l 



Ws{x)\'^ (6) 



where q includes only lowest-band states since A^ < 1^. 
Hence, the BHM is valid for the whole range of inter- 
actions. We can understand the final profile also in a 
hole-excitation picture where starting from the the MI 
in the commensurate case N = W we annihilate W — N 
particles ^ = X]a=i cia\MlN^w)- While in the case of 
v = 1 the addition of A^ = 1^ coefficients in Eq. 6 leads 
to equal site occupation, here we can have imbalances 
and oscillations of the density depending on how many 
orbitals contribute according to the numbers W and N 
(Fig. [TOTa) g = 3.0, Fig. [TOTb) inset). Triggered by the 
interaction, the 'transfer' of particles from the middle 
of the potential to the outer positions passes through 
the intermediate wells which gain and lose population 
(Fig. [TOT b) inset 5 = 2). The number fluctuations in 
Fig. [To] (b) saturate to a rather high value because of 
the incommensurate filling u < 1 which allows only delo- 
calised phases. The fluctuations are greater in the wells 
with less population corresponding to 'holes' (Fig.[TQl(b) 
compare s = 2, 3 with s = 1,4). 

A main difference between commensurate and incom- 
mensurate filling is that in the latter case the 'coherence', 
or better the off-diagonal part of pi(x, x'), cannot vanish 
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Figure 10: (Colour online) (a) p{x) for 5 particles and 7 wells. 
Shown are 4 different values of g: non interacting (g = 0-0), 
weakly interacting (g = 0.5,^ = 1.0) and fermionized limit 
(g — 3.0). (b) Particle number fluctuations as g increases for 
the sites s — 1,2,3,4. (inset) On-site population. 



completely since the particles remain in fact delocalised. 
The remaining one-body correlations in the fermioniza- 
tion limit (Fig. [11] (a)) are concentrated mostly close to 
the diagonal, i.e., between neighbouring lattice sites. It 
is interesting though that this short range coherence is 
not equally distributed on all close to the diagonal spots. 
For example, in the case of 5 particles in 7 wells, there 
is a quite well localised particle in the central site as the 
vanishing off-diagonal terms in the center of the figure 
indicate. The distribution of short-range correlations is 
again related to the commensur ability of the setup, and 
can be understood in the following way: divide the real 
space of the multi-well potential into N equal intervals 
and put each particle in the middle of one interval; then 
those which lie closer to the middle of a site are also 
better localised than those who lie close to the inter-well 
barriers affecting the one-body correlations accordingly, 
i.e., when there is more localisation the correlations die 
out. 

The attempted localisation of the particles as the re- 
pulsion increases distorts to some extent the interference 
pattern in the momentum distribution (Fig. [11] (c)). In 
particular the central peak is lowered due to the partial 
loss of coherence but the peaked structure is not fully 
smeared out. We do not observe an increase of the cen- 
tral peak of the momentum for low interactions as in the 
case oi V = 1, even though we use a rather shallow lat- 
tice Vq = 10.0 = 4.1^;^^ with d = 1.42). Although the 
first natural orbital is broadened and dominant [Fig. [11] 
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Figure 11: (Colour online) (a) pi{x, x') for 7 wells and 5 parti- 
cles in the fermionization limit, (b) Population of the natural 
orbitals as a function of g. (c) Momentum distribution for 5 
different values of g: non interacting {g — 0.0), weakly inter- 
acting [g — 0.1), fermionization limit {g — 3.0, 20.0). 



(b)], the density of 'condensed' particles ^^^ is not as 
high in this case of incommensurate filling i^ < 1 as for 
the case v = 1 and this explains the instant lowering of 
the A: = peak. As expected W natural orbitals con- 
tribute substantially (see Fig. [TT] (b)) as for i/ = 1, but 
the contribution of each orbital here, differs throughout 
the fermionization crossover. 



B. Incommensurate Filling v > 1 

The incommensurate v > 1 case combines localisation 
delocalisation as well as on-site interaction effects. It is 
instructive to keep in mind the properties of the corre- 
sponding v^ = Nrn^w f^ii^j^g case, which refers only to 
the 'extra' NmodW particles, thereby pointing out sim- 
ilarities and differences which allow to identify the effect 
of the 'background'. For all the following cases the length 
unit is L/9 and d = 2.2 (four wells). The lattice depth 
is chosen Vb = 20.0 = 19.6 Er such that the particles are 
confined energetically below the continuum. 



1. One extra particle on a unit filling background 

We begin with the simplest case of one extra particle 
added to the unit filling (here 5 particles in 4 wells). For 
weak interactions the ^ = nonuniform occupation of 
sites tends to become uniform [Fig. [T2ja) g = 0.05]. For 
a slightly higher interaction {g = 1.0) strength though. 



there is an interesting revival of the tendency to predom- 
inantly occupy the middle wells. To understand this, 
we need to go beyond the lowest band and BHM anal- 
ysis (since u > 1) and consider contributions of higher 
bands. The higher band states possess a similar distri- 
bution with respect to the different wells as the lowest 
band (see Eq. 2 ). Thus the energetically lowest level 
of the excited band, which has a dominant population 
in the middle (see eg. Fig. [1] for the triple well), when 
contributing, results in repopulation of the center. Ad- 
ditionally as the interaction increases, the total energy 
of the particles becomes higher and thus they approach 
energetically the top of the barriers of the potential close 
to the continuum (see also the case u = 2). This en- 
forces the hopping term since higher bands with larger 
coupling J^ > J^ contribute (compare intra-band split- 
tings in Fig. [1]) and as a consequence the kinetic energy 
term redirects the particles to the center. This intuitive 
picture of contributions from higher energy states that 
we have drawn here for the repopulation of the middle, 
is also consistent with the treatment of the Tonks limit 
via the Bose-Fermi map. According to the theorem, the 
extra particle lies exactly on the energetically lowest one 
particle level of the excited band, which possesses also 
higher contributions in the middle, while the other par- 
ticles completely occupy the lowest band states forming 
a MI background of one particle localised per well. This 
results in a slight broadening especially of the central 
peaks of the one-body density (Fig. [I2fa) g = 20.0), a 
standard on-site interaction effect for v > 1^ which we 
encountered also in the previous section for v = 2. Let 
us note that the situation is quite different for a suf- 
ficiently shallow lattice not considered here: the extra 
particle would go closer or even above the barriers and 
thus would fill the inter- well space, distributing smoothly 
over the potential and resulting in strongly reduced on- 
site effects [33 • A reference model to understand incom- 
mensurate filling i^ > 1, in qualitative agreement with 
our results, was given in [51|: the particles occupy dif- 
ferent horizontal 'layers', each one on top of the other, 
all having commensurate MI states and only the highest 
one being incommensurate (with NvnodW particles) and 
delocalised. 

In the evolution of the populations with increasing g 
we can verify the density variation in the center wells (see 
Fig. [T2r b) inset). In the strongly interacting regime, the 
populations remain quite similar due to the background 
of localised particles and they differ only because of the 
non-uniform distribution of the extra particle in the first 
level of the excited band. The particle-number fluctu- 
ations (Fig. [T2r b)) remain quite large, since the extra 
delocalised particle does not allow for a perfect insulator 
phase. Nevertheless they are substantially diminished 
compared to the corresponding v'^ =1/4 single particle 
case because of the localisation of the background. 

The non-local properties confirm the incomplete lo- 
calisation in this case. In the one body density matrix 
pi{x^x') the remaining coherence in the strongly inter- 
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Figure 12: (Colour online) (a) p{x) for 5 particles and 4 wells. 
Shown are 4 different values of g: non interacting (g — 0.0), 
weakly interacting [g — 0.05,1.0,5.0), fermionization limit 
{g — 20.0). (b) Particle number fluctuations as g increases 
for the sites s — 1,2. (inset) On-site populations. 
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Figure 13: (Colour online) (a) One-particle density matrix 
p(x, x') for 5 particles in 4 wells in the non-interacting g — 0.0 
and (b) in the fermionization limit g — 20.0. (c) Mo- 
mentum distribution. Shown are 5 different values of g: 
(g = 0.0,0.05,1.0,5.0,20.0). (d) Population of the natural 
orbit als as a function of the interaction strength. 



acting limit (Fig. [13] (b)), is concentrated close to the di- 
agonal. Due to the localised background the long range 
off-diagonal terms almost disappear, compared to the 
^ = 1/4 case of a single particle (Fig. [13] (a)). The off- 



u 



diagonal humps, which are mainly visible in the center, 
reflect a widened pattern (Fig. [13] (b)). In the momentum 
distribution (Fig. [TST c)) the central peak and the Bragg 
peaks are from very weak interactions {g = 0.05, 1.0) low- 
ered because of the large depth of the potential used in 
this case {Vq = 20.0 = 19.6^i^). Small peaks still persist, 
revealing incomplete localisation due to the extra particle 
(Fig.[l3l^c) for ^ = 20.0). 



2. A repulsive pair of particles on a localised background 

Choosing the number of particles and wells at will, one 
encounters many different incommensurate cases and cor- 
responding effects. We focus here on the case of incom- 
mensurate filling with two extra particles, which exhibits 
a distinct behaviour: the extra pair of particles feels the 
interaction with the background particles, additionally 
to the intra-pair repulsive forces. From the above discus- 
sion (sec IV. B.l), we recall tendencies to first equalise 
and then repopulate the center wells with increasing g. 
Indeed this happens also in the low-interaction regime 
for the case of 6 particles in 4 wells examined here (see 
Fig. [T4lf a) g = 0.1, 1.0, Fig. [T4lf c)). It even results in a 
broadening of the central peaks and wiggles in the one- 



body density for higher interaction strengths (Fig. [MT a) 
g = 10.0,30.0). To realise the peculiarity of this ef- 
fect one has to consider the fermionization limit which 
exhibits only a slight broadening but no wiggles at all; 
rather the particles are distributed uniformly (Fig. [T4lf a) 
and Fig. [TJ] (b) for g = 100.0) in accordance to the 
Bose-Fermi map which predicts equal contribution from 
the two first levels of the excited band. However, for 
strong but finite interactions the extra particles concen- 
trate more in the center resulting in an effectively higher 
local filling 1^3=2,3 ^ 2; only for ^ > 30 the repulsion is 
strong enough to drive the bosons to the outer wells. 

In Fig. [T4t^ c) we observe that for weak to intermediate 
interactions {g ~ 0.05—0.8), the populations and the fluc- 
tuations are almost constant. Beyond this regime (close 
to ^ = 0.8), the population of the center wells {s = 2) 
increases and approaches the value for two particles in 
these sites 1^3=2,3 ^ 2, while the fluctuations are again 
strongly reduced. This can be understood as a localisa- 
tion behaviour ('Mott-like phase') in a central 'domain' of 
the potential in analogy with similar situations appearing 
for an optical lattice with a superimposed harmonic con- 
finement 0, S [IqI, [i3-[il Ell-Is^ • in the latter case, the 
harmonic potential increases the on-site energy as we go 
to the outer wells and thus sets an 'energetic obstacle' for 
the particles to occupy them. Therefore, they prefer to 
localise in the center and form Mott domains (or shells), 
possibly surrounded by a superfluid layer of delocalised 
particles. As an illustration of this effect in our few-body 
setup, we present the case of 5 particles in a ID lattice 
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with a superimposed harmonic trap (Fig. [MT a) inset). 
For strong interactions {g = 30.0 here) we have exactly 
two particles in each middle well and one particle divided 
into the two outer wells surrounding the 'Mott-shelP. We 
underline that our numerically exact method which goes 
beyond BHM brings a new light on the strongly interact- 
ing regime including the on-site interaction effects; here 
we point out the formation of wiggles in the one-body 
density inside the 'Mott-shell' of two particles per site. 
For our initial setup though there is no harmonic con- 
finement and thus no 'energetic obstacles' between the 
wells; nevertheless the hard- wall boundary conditions on 
the edges, make it preferable for the particles to be in 
the center. We can thus comment that the finite size 
itself makes it possible for qualitatively different spatial 
regions to occur in incommensurate one-dimensional lat- 
tices. For stronger interactions, the fluctuations increase 
again (Fig. [ll](b)) as the particles occupy higher energy 
levels and delocalise further. The fermionization comes 



with equal population -^ = 1.5 particles per well (Fig.fH 



(b) inset). 

Similar phenomena of course can happen for other 
cases of incommensurate filling u > 1. Let us stress that 
the exact number of wells and particles is important for 
the possible effects. For example, in the case of 5 particles 
in 3 wells (Fig. [14] (d)), the repulsively interacting extra 
pair of particles is prohibited from occupying the middle 
well, since it would require a very undesirable triple oc- 
cupation of a single site with strongly repulsive particles. 
Thus the revival of occupation in the center is avoided, 
and each of the two extra particles is mostly located in 
one outer well, indicated by the corresponding wiggles of 
the one-body density. 

Let us proceed by analysing the case of 6 particles 
in 4 wells from a two-body perspective. For weak in- 
teractions, the off-diagonal peaks of ^2(^1,^2) become 
more pronounced than the diagonal contribution (Fig. [15] 
g = 0.1). For stronger interactions {g = 10.0) we ob- 
serve reconcentration in the center and the correlation 
hole is starting to be formed. The off-diagonal humps 
start to broaden in the middle and acquire fragmented 
patterns while the correlation hole becomes more evident 
{g = 30.0). In the fermionization limit {g = 100.0) the 
distribution on the diagonal becomes equal for all sites, 
and the off-diagonal contributions lose their fragmented 
pattern. As a last comment for the comparison of the 
crossover seen in the one- and two-body density, we ob- 
serve that the formation of correlation hole in the diag- 
onal of p2{xi^X2) (Fig. [I5] ^ = 10.0) corresponds to a 
widening of the corresponding peaks in p{x) (Fig. [HT a) 
g = 10.0), while the appearance of wiggles (Fig. [MT a) 
g = 30.0) is connected with the fragmented patterns 
in the corresponding off-diagonal humps of P2(^i,^2) 
(Fig. [15]^ = 30.0). 

An interesting effect occurs in this case for the one- 
body density matrix pi{xjx'). The remaining coherence 
in the strongly interacting limit is only concentrated in 
the left and the right part of the space (diagonal squares 
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Figure 14: (Colour online) (a) One body-density for 6 parti- 
cles and 4 wells. Shown are 6 different values of g: non inter- 
acting (g — 0.0), weakly interacting (^ = 0.1, 1.0), strongly in- 
teracting (g — 10.0, 30.0) and fermionization limit [g — 10.0). 
(inset) The case of a superimposed harmonic trap with strong 
interactions g — 30.0 for 5 particles, (b) Particle number 
fluctuations as a function of g for s — 1,2. (inset) On- 
site populations, (c) Same plot for the weak interaction 
regime, (d) One-body density for 5 particles in 3 wells for 
^ = 0.0,0.05,5.0,20.0. 



in Fig. [16] (a), (b)). This indicates that the two extra par- 
ticles, tend to localise with respect to each other since the 
off-diagonal squares of this plot are completely depleted; 
the left is completely uncorrelated with the right part of 
the space. Apparent is the formation of wiggles in the 
center of the diagonal for g = 30.0 (Fig.[T6]^a)) along with 
the formation of fragmented patterns in the remaining 
off-diagonal contribution. For the fermionization limit 
(Fig. [TCT b) g — 100.0) these wiggles are smeared out to an 
equally broadened profile, while there is a very strongly 
fragmented pattern in the populated off-diagonal spots. 
It is worth comparing this to the fermionization limit 
with the corresponding v'^ = 2/4 case, 2 particles in 4 
wells (Fig. [16] (c)). The coherence between left and right 



14 




(c) 



-4-3-2-1 1 2 3 4 



(a) 


•^^ 


4 
3 


• «« 


2 

1 

-1 
-2 


1 A ^# i 


-3 
-4 


•• • 



-4-3-2-1 1 2 3 4 




-4-3-2-1 1 2 3 4 

(d) x^ 



Figure 15: (Colour online) Two body-density p2ixi,X2) for 6 
particles and 4 wells for (a) ^ = 0.1, (b) g = 10.0, (c) g — 30.0, 
(d) g = 100.0. 
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Figure 16: (Colour online) One-body density matrix p{x^x') 
for 6 particles in 4 wells (a) g — 30.0, (b) g — 100. and (c) 
2 particles in 4 wells fermionization limit g — 20.0. (d) Mo- 
mentum distribution for 6 particles in 4 wells. Shown are 6 
different values of ^: non interacting [g — 0.0), weak interact- 
ing {g = 0.005,0.1,2.0), strong interactions {g = 10.0,30.0) 
and fermionization limit [g — 100.0) which is compared with 
the fermionization limit of 2 particles in 4 wells. 



part of the potential is present in the latter case (off- 
diagonal squares retain contribution), which means that 
one particle mainly localised in the left part does pene- 
trate into the right part. Contrarily, the background of 
localised particles in the case of 6 particles in 4 wells, pre- 
vents each of the excited extra particles from intruding 
into the side where the other one sits. 

The localisation of the extra particles in the two 
middle wells for strong interactions leads to an almost 
smoothened (MI) profile of the momentum distribution 
(Fig. fTGT d) g = 10.0). As we go to the fermioniza- 
tion limit {g = 100.0), the central peak becomes a little 
higher, and the profile is somewhat distorted again, as 
the coherence slightly increases because of the delocali- 
sation of the extra particles. Of course the corresponding 
v^ = 2/4 case (Fig. [TCT d) 2 particles in 4 wells) is much 
less smoothened, since the delocalisation is not hindered 
by any localised background. 



V. CONCLUSIONS AND OUTLOOK 

We have performed a numerically exact investigation 
of few-boson systems in finite one-dimensional lattices 
for varying strength of the repulsive interactions. In the 
unit filling case, we have shown the evolution from a 
non-interacting state with a maximum of the density in 
the center to an equal site distribution and a simulta- 
neous decrease of the fluctuations; these effects are in 
accordance with the predictions of the Bose Hubbard 
model, following the Superfluid - Mott insulator tran- 



sition from a delocalised to a localised state. Beyond 
that, for higher commensurate filling, on-site interaction 
effects like broadening and wiggles in the one-body den- 
sity as well as correlation hole and fragmented patterns 
in the two body density occur especially in the strongly 
repulsive limit approaching fermionization. The coher- 
ence loss due to the interaction is reflected in the re- 
duced off-diagonal contribution in the one-body density 
matrix and the smoothening of the peaked pattern in 
the momentum distribution. The high-momentum tails 
show a difference between the Mott-insulator within the 
Bose-Hubbard model with unperturbed on-site functions 
and the Tonks limit with two fermionized particles per 
site. Interestingly, there is a slight increase of long-range 
correlations on the onset of interactions with the parti- 
cles remaining mostly in one orbital which broadens. The 
fragmentation into different natural orbitals for strong in- 
teractions reflects the band structure of the lattice. The 
effect of a deeper lattice on enhancing and accelerating 
these processes was pointed out. For incommensurate fill- 
ing we have shown that the density distributes inhomoge- 
neously depending on the number of particles and wells, 
and partially delocalised incoherent states are formed. In 
particular for filling factor greater than one, the degree of 
localisation of these extra particles depends on the inter- 
action strength, the presence of the background particles 
as well as the specific conditions of the setup. For finite 
systems, it can vary locally, leading to spatial variations 
for the observables, in analogy with the insulating and 
superfluid domains which appear in the case of an addi- 
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tional harmonic confinement. Our study suggests many 
promising routes for further investigations, like spatially 
inhomogeneous lattices, particularly the case of an ex- 
ternal harmonic trap and disorder with arbitrary energy 
offsets between the wells. 
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Dirac-Frenkel variational principle, one can derive equa- 
tions of motion for both Aj, ^j. This conceptual compli- 
cation offers an enormous advantage: the basis {<l>j(-,t)} 
is variationally optimal at each time t, allowing us to 
keep it fairly small. The permutation symmetry can be 
enforced by symmetrising the coefficients Aj, but the 
ground state is automatically bosonic. 

The Heidelberg MCTDH package [59], incorporates 
the so-called relaxation method which provides a way to 
obtain the lowest eigenstates of the system by propa- 
gating some wave function ^o by the non-unitary e~^^ 
{propagation in imaginary time.) As r ^ oo, this auto- 
matically damps out any contribution but that stemming 
from the true ground state |0), 



Appendix A: Computational Method 

The goal is to investigate the ground state properties of 
finite bosonic systems introduced in Sec. H for the com- 
plete range of interaction strengths in a numerically exact 
way. Our approach relies on the Multi-Configurational 
Time-Dependent Hartree (MCTDH) method [45 , IE, [58J , 
primarily a wave-packet dynamics tool known for its 
outstanding efficiency in high-dimensional applications. 
The underlying idea of MCTDH is to solve the time- 
dependent Schrodinger equation 



i^ = H^ 

^(Q,o) = ^o(g) 



(Al) 



as an initial- value problem by expansion in terms of di- 
rect (or Hartree)products $j: 



^(Q,t) = ^Aj(t)$j(Q,t) 



J 



ni 



n 



/ 






/^=1 



using a convenient multi-index notation for the configu- 
rations, J = (ji . . . j/), where f = N denotes the number 
of degrees of freedom and Q = (xi, . . . ,Xf)^ . The single- 
particle functions Lp y' are in turn represented in a fixed, 
primitive basis implemented on a grid. For indistinguish- 
able particles as in our case, the single-particle functions 
for each degree of freedom /^ = 1 , . . . , A/^ are of course 
identical in both type and number {^j^^ with j^ < n). 

In the above expansion, both the coefficients Aj and 
the Hartree products ^j are time-dependent. Using the 



-Ht 



*o = E 



-Ejr 



^)(^l^o) 



J 



In practice, one relies on a more sophisticated scheme 
termed improved relaxation [60]. Here {^\H — E\^) is 
minimised with respect to both the coefficients Aj and 
the configurations <l>j. The equations of motion are 
solved iteratively, first for Aj(t) (by diagonalisation of 
{{^j\H\^k)) with fixed ^j) and then propagating ^j 
in imaginary time over a short period. The cycle will 
then be repeated. 

As it stands, the effort of this method scales exponen- 
tially with the number of degrees of freedom, n^ . Just 
as an illustration, using 15 orbitals and A^ = 5 requires 
7.6 • 10^ configurations J. This restricts our analysis in 
the current setup to about TV = 0(10), depending on 
how decisive correlation effects are. If these are indeed 
essential, then it turns out that at least n = N orbitals 
are needed for qualitative convergence alone, while the 
true behaviour may necessitate about 15 for A^ = 5. By 
contrast, the dependence on the primitive basis, and thus 
on the grid points, is not as severe. In our case, the grid 
spacing should of course be small enough to sample the 
interaction potential, and we consider a basis set of sinu- 
soidal functions which guarantee the hard- wall boundary 
condition (zero value of the wave function on the first and 
the last grid point). Of interest from a methodological 
point of view (but also from a conceptual one) is that a 
shallower lattice enhances the convergence of this com- 
putational method. For a deep lattice the localisation of 
the particles in individual wells is pronounced and thus 
a higher number of delocalised single-particle functions 
are necessary to express them. 
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